compare.size.frq <- function(frq1,frq2,fishery=5, wt=T, prefx = "_",doyears, fdesc="",summary=TRUE)
{
##==================================================================================
## takes two frq files and compares the length and weight data on an annual basis for
##  the range of years -- just does the last 42 years at the moment
##
## modified by SDH July 2010 to generalize and allow for no wt data
## NMD 30 May 2011 - modified to:
## - take account of frq file version, 
## - plots proportions not numbers,
## - includes an argument "prefx_" to give output plot file a unique name
## - returns a list object for each frq file with summaries of quantities
## - deals with frq files having different numbers of fisheries
## List structure: fish=fishery,len=out.len,wei=out.wei,catch=ann.catch,effort=ann.effort,
##               : cpue=ann.cpue,out.len and out.wei: year, p25, p50, p75, sample_size
##
## still needs work to generalize the weight data section, handle different frq versions
##==================================================================================

  # Get version number of frq file to determine first column for checking data
  vsn1 <- frq1$struct$te
  vsn2 <- frq2$struct$te
  if(vsn1 != vsn2){ 
    print("Freq files input are different versions!")
    break()
  }
  col1 <- ifelse(vsn1 < 6, 7, 8)          # First column for length frequencies
  
  # Check numbers of fisheries in each frq object
  pltfrq1flg <- 0
  pltfrq2flg <- 0
  fadd <- 0
  frq1n <- frq1$struct$nf
  frq2n <- frq2$struct$nf
  if(frq1n != frq2n){
    nfdif <- abs(frq1n - frq2n)
    if(frq1n > frq2n){    # add dummy fisheries to frq2
      fadd <- seq(from = frq2n+1, to = frq1n, by = 1)
      tmp <- frq1$mat[frq1$mat[,4] %in% fadd,]        # copy missing fisheries from frq1
      tmp[,c(col1,(col1+1))] <- -1
      tmp[,c((col1+2):dim(frq2$mat)[2])] <- 0
      frq2$mat <- rbind(frq2$mat,tmp)
      pltfrq2flg <- 1
    } else {    # add dummy fisheries to frq1
      fadd <- seq(from = frq1n+1, to = frq2n, by = 1)
      tmp <- frq2$mat[frq2$mat[,4] %in% fadd,]        # copy missing fisheries from frq2
      tmp[,c(col1,(col1+1))] <- -1
      tmp[,c((col1+2):dim(frq1$mat)[2])] <- 0
      frq1$mat <- rbind(frq1$mat,tmp)
      pltfrq1flg <- 1
    }
  }
  ## get data
  #frq1 <- yyy
  data1 <- frq1$mat[frq1$mat[,4]==fishery,]
  # frq2 <- xxx
  data2 <- frq2$mat[frq2$mat[,4]==fishery,]
  
  lfint <- frq1$dl$lfint
  wfint <- frq1$dl$wfint
  lffirst <- frq1$dl$lffirst
  wffirst <- frq1$dl$wffirst
  lfwidth <- frq1$dl$lfwidth
  wfwidth <- frq1$dl$wfwidth
  
  lens <- seq(from=lffirst,length=lfint,by=lfwidth)
  weis <- seq(from=wffirst,length=wfint,by=wfwidth)
  lcol <- lfint+2
  wcol <- wfint+2
  
  # first create objects for length and weight for each dataset - missing values are zeros
  len.1 <- matrix(0,nrow=nrow(data1),ncol=lcol)
  len.1[,1:2] <- data1[,1:2]
  wei.1<- matrix(0,nrow=nrow(data1),ncol=wcol)
  wei.1[,1:2] <- data1[,1:2]
  # first data set
  for(i in 1:nrow(data1))
  {
    #Length data first
    if(data1[i,col1]==-1)
    {
    # do nothing as already zeros
      if(wt==F)
      {
      # no weight data either
      }
      else if (data1[i,(col1+1)]==-1)
      {
      }
      else
      {
      wei.1[i,3:wcol] <- data1[i,(col1+1):(col1+wcol-2)]
      }
    }
    else
    {
    len.1[i,3:(2+lfint)] <- data1[i,col1:(col1-1+lfint)]
      if(wt==F)   # col1 -1 gets to start of LF,
                                                  # + length of LF intervals (lfint)
                                                  # +1 to position of WF pointer
      {
      # no weight data
      }
      else if(data1[i,(col1-1+lfint+1)]==-1)
      {
      }
      else
      {
      wei.1[i,3:wcol] <- data1[i,(col1-1+lfint+1):(col1-1+lfint+wfint)]  # col1 -1 gets to start of LF,
                                                                         # + length of LF intervals (lfint)
                                                                         # +1 to position of WF pointer
      }
    }
  }
  # second dataset
  # first create an objects for length and weight for each dataset - missing values are zeros
  len.2 <- matrix(0,nrow=nrow(data2),ncol=2+lfint)
  len.2[,1:2] <- data2[,1:2]
  wei.2<- matrix(0,nrow=nrow(data2),ncol=wcol)
  wei.2[,1:2] <- data2[,1:2]
        for(i in 1:nrow(data2))
        {
            #Length data first
            if(data2[i,col1]==-1)
            {
            # do nothing as already zeros
                if(wt==F)
                {
                # no weight data either
                }
                else if (data2[i,(col1+1)]==-1)
                {
                # no weight data either
                }
                else
                {
                wei.2[i,3:wcol] <- data2[i,(col1+1):(col1+wcol-2)]
                }
            }
            else
            {
            len.2[i,3:(2+lfint)] <- data2[i,col1:(col1-1+lfint)]
                if(wt==F)             # col1 -1 gets to start of LF, + length of LF intervals (lfint) +1 to position of WF pointer
                {
                # no weight data
                }
                else if(data2[i,(col1-1+lfint+1)]==-1)
                {
                }
                else
                {
                wei.2[i,3:wcol] <- data2[i,(col1-1+lfint+1):(col1-1+lfint+wfint)]    # col1 -1 gets to start of LF, + length of LF intervals (lfint) +1 to position of WF pointer
                }
            }
        }
  
  
  
  # Length plot
  if(!missing(doyears)){
   years <- doyears
  } else {
    years <- unique(c(data1[,1],data2[,1]))
  }
      # if more than 42 years then just do the last 42
      if(length(years)>42){
       p.years <- rev(rev(years)[1:42])
      } else {
      p.years <- years
      }
  
      #browser()
  
  ##x11(height=17,width=28)
  ncolpl <- ceiling(sqrt(length(p.years)))
  nrowpl <- ceiling(length(p.years)/ncolpl)
  par(mfcol=c(nrowpl,ncolpl),mar=c(1,1,1,1),oma=c(0.5,0.5,2,0.5))
  #par(mfcol=c(6,7),mar=c(1,1,1,1),oma=c(0.5,0.5,2,0.5))
      for(i in 1:length(p.years))
      {
        nlen1 <- 0        # Initialise total sample sizes
        nlen2 <- 0
        #print(i)
        tmp.1 <- len.1[len.1[,1]==p.years[i],]
        #browser()
        if(is.null(nrow(tmp.1))){
          len.sum.1 <- tmp.1[3:(2+lfint)]
          len.sum.1[is.na(len.sum.1)] <- 0 
          if(sum(len.sum.1) > 0.0001){
            nlen1 <- sum(len.sum.1)
            len.sum.1 <- len.sum.1 / sum(len.sum.1)
          }
          len.sum.1[len.sum.1 == Inf | len.sum.1 == -Inf] <- 0
        } else {
          len.sum.1 <- apply(tmp.1[,3:(2+lfint)],2,sum)
          len.sum.1[is.na(len.sum.1)] <- 0 
          if(sum(len.sum.1) > 0.0001){
            nlen1 <- sum(len.sum.1)
            len.sum.1 <- len.sum.1 / sum(len.sum.1)
          }  
          len.sum.1[len.sum.1 == Inf | len.sum.1 == -Inf] <- 0
        }
    
        tmp.2 <- len.2[len.2[,1]==p.years[i],]
        if(is.null(nrow(tmp.2))){
          len.sum.2 <- tmp.2[3:(2+lfint)]
          len.sum.2[is.na(len.sum.2)] <- 0 
          if(sum(len.sum.2) > 0.0001){
            nlen2 <- sum(len.sum.2)
            len.sum.2 <- len.sum.2 / sum(len.sum.2)
          }
          len.sum.2[len.sum.2 == Inf | len.sum.2 == -Inf] <- 0
        } else {
          len.sum.2 <- apply(tmp.2[,3:(2+lfint)],2,sum)
          len.sum.2[is.na(len.sum.2)] <- 0 
          if(sum(len.sum.2) > 0.0001){
            nlen2 <- sum(len.sum.2)
            len.sum.2 <- len.sum.2 / sum(len.sum.2)  
          }
          len.sum.2[len.sum.2 == Inf | len.sum.2 == -Inf] <- 0
        }
         ymax <- max(c(len.sum.1,len.sum.2))
        plot(lens,len.sum.1,ylim=range(c(0,ymax)),type="n",xlab="",ylab="",axes=F)
        lines(lens,len.sum.1,lwd=2,col=1,lty=1)
        lines(lens,len.sum.2,lwd=2,col=2,lty=2)
        box()
        mtext(side=3,outer=F,text=p.years[i],cex=0.7)
        legend(10,ymax,legend=c("Frq1","Frq2"),lty=c(1,1),lwd=c(2,2),col=c(1,2))
        text(160, (0.9*ymax) , paste("n_Frq1= ",nlen1,sep=""),  adj = NULL, pos = NULL,
             offset = 0.5, vfont = NULL, cex = 1, col = 1, font = NULL)
        text(160, (0.8*ymax) , paste("n_Frq2= ",nlen2,sep=""),  adj = NULL, pos = NULL,
             offset = 0.5, vfont = NULL, cex = 1, col = 2, font = NULL)
      }
      if(!fishery %in% fadd){          # frq files have same number of fisheries
        ##browser()
        mtext(side=3,outer=T,text=paste("length-",fishery," ",fdesc[fishery,2],sep=""))
      } else {            # add note to plot for missing fisheries
        if(pltfrq1flg){
          mtext(side=3,outer=T,text=paste("length",fishery," - fishery is missing in frq1",sep="-"))
        }
        if(pltfrq2flg){
          mtext(side=3,outer=T,text=paste("length",fishery," - fishery is missing in frq2",sep="-"))
        }
      }
      if(Sys.info()[1] != "Linux") {
        savePlot(filename=paste(prefx,"_length",formatC(fishery,width=2,flag="0"),".png",sep="_"),type="png")
      } else {
        pngfile(paste(prefx,"_length",fishery,sep="-"))
      }
          
  
  #  Weight frequency plots
  if(wt) {
    par(mfcol=c(ncolpl,nrowpl),mar=c(1,1,1,1))
    #par(mfcol=c(6,7),mar=c(1,1,1,1))
      for(i in 1:length(p.years))
      {
        nwt1 <- 0        # Initialise total sample sizes
        nwt2 <- 0
        tmp.1 <- wei.1[wei.1[,1]==p.years[i],]
        if(is.null(nrow(tmp.1))){
          wei.sum.1 <- tmp.1[3:wcol]
          wei.sum.1[is.na(wei.sum.1)] <- 0 
          if(sum(wei.sum.1) > 0.0001){
            nwt1 <- sum(wei.sum.1)
            wei.sum.1 <- wei.sum.1 / sum(wei.sum.1)
          }
          wei.sum.1[wei.sum.1 == Inf | wei.sum.1 == -Inf] <- 0
        } else {
          wei.sum.1 <- apply(tmp.1[,3:wcol],2,sum)
          wei.sum.1[is.na(wei.sum.1)] <- 0 
          if(sum(wei.sum.1) > 0.0001){
            nwt1 <- sum(wei.sum.1)
            wei.sum.1 <- wei.sum.1 / sum(wei.sum.1)
          }
          wei.sum.1[wei.sum.1 == Inf | wei.sum.1 == -Inf] <- 0
        }
      
        tmp.2 <- wei.2[wei.2[,1]==p.years[i],]
        if(is.null(nrow(tmp.2))){
          wei.sum.2 <- tmp.2[3:wcol]
          wei.sum.2[is.na(wei.sum.2)] <- 0 
          if(sum(wei.sum.2) > 0.0001){
            nwt2 <- sum(wei.sum.2)
            wei.sum.2 <- wei.sum.2 / sum(wei.sum.2)
          }
          wei.sum.2[wei.sum.2 == Inf | wei.sum.2 == -Inf] <- 0
        } else {
          wei.sum.2 <- apply(tmp.2[,3:wcol],2,sum)
          wei.sum.2[is.na(wei.sum.2)] <- 0 
          if(sum(wei.sum.2) > 0.0001){
            nwt2 <- sum(wei.sum.2)
            wei.sum.2 <- wei.sum.2 / sum(wei.sum.2)
          }
          wei.sum.2[wei.sum.2 == Inf | wei.sum.2 == -Inf] <- 0
        }
    
        ymax <- max(c(wei.sum.1,wei.sum.2))
    
        plot(weis,wei.sum.1,ylim=range(c(0,ymax)),type="n",xlab="",ylab="",axes=F)
        lines(weis,wei.sum.1,lwd=2,col=1,lty=1)
        lines(weis,wei.sum.2,lwd=2,col=2,lty=2)
        box()
        mtext(side=3,outer=F,text=p.years[i],cex=0.7)
        legend(10,ymax,legend=c("Frq1","Frq2"),lty=c(1,1),lwd=c(2,2),col=c(1,2),bty="n")
        text(160, (0.9*ymax) , paste("n_Frq1= ",nwt1,sep=""),  adj = NULL, pos = NULL,
             offset = 0.5, vfont = NULL, cex = 1, col = 1, font = NULL)
        text(160, (0.8*ymax) , paste("n_Frq2= ",nwt2,sep=""),  adj = NULL, pos = NULL,
             offset = 0.5, vfont = NULL, cex = 1, col = 2, font = NULL)
      }
      if(!fishery %in% fadd){          # frq files have same number of fisheries
        mtext(side=3,outer=T,text=paste("weight-",fishery," ",fdesc[fishery,2],sep=""))
      } else {            # add note to plot for missing fisheries
        if(pltfrq1flg){
          mtext(side=3,outer=T,text=paste("weight",fishery," - fishery is missing in frq1",sep="-"))
        }
        if(pltfrq2flg){
          mtext(side=3,outer=T,text=paste("weight",fishery," - fishery is missing in frq2",sep="-"))
        }
      }
   
      if(Sys.info()[1] != "Linux") {
        savePlot(filename=paste(prefx,"_weight",fishery,sep="-"),type="png")
      } else {
        jpegfile(paste(prefx,"_weight",fishery,sep="-"))
      }
  }
  # Now create summary list objects for quantities in frq files
  if(summary) {
    compr_frq <- list()
    compr_frq[[1]] <- list()
    compr_frq[[2]] <- list()
    lstnames <- c("frq1_lst","frq2_lst")
    names(compr_frq) <- lstnames
    
    # Identify frq objects
    lenfrq <- c("len.1","len.2")
    weifrq <- c("wei.1","wei.2")
    datfrq <- c("data1","data2")
    
    for(i in 1:length(lstnames)){
    # aggregate to the year
      ann.len <- aggregate(get(lenfrq[i])[,3:(2+lfint)],list(get(lenfrq[i])[,1]),sum)
      if(wt) ann.wei <- aggregate(get(weifrq[i])[,3:wcol],list(get(weifrq[i])[,1]),sum)
      ann.catch <- aggregate(get(datfrq[i])[,5],list(get(datfrq[i])[,1]),sum)
      ann.effort <- aggregate(get(datfrq[i])[,6],list(get(datfrq[i])[,1]),sum)
      ann.cpue <- ann.effort
      ann.cpue[,2] <- ann.catch[,2]/ann.effort[,2]
      
      # setup data.frame to take quantiles for length and weight
      out.len <- data.frame(year=ann.len[,1],
                            p25=rep(NA,length=nrow(ann.len)),
                            p50=rep(NA,length=nrow(ann.len)),
                            p75=rep(NA,length=nrow(ann.len)), samp=NA)
      if(wt) out.wei <- out.len
      
      # get the quantiles of length and weight for each year
      for(j in 1:nrow(ann.len))
      {
        if(rowSums(ann.len[j,2:lfint+1])!=0)
        {
          out.len[j,2:4] <- quantile(rep(lens,ann.len[j,1:lfint+1]),probs=c(0.25,0.5,0.75))
          out.len[j,5] <- rowSums(ann.len[j,1:lfint+1])
        }
      
        if(wt) if(rowSums(ann.wei[j,2:lfint+1])!=0)
        {
          out.wei[j,2:4] <- quantile(rep(weis,ann.wei[j,2:(wcol-1)]),probs=c(0.25,0.5,0.75))
          out.wei[j,5] <- rowSums(ann.wei[j,2:lfint+1])
        }
      } 
      if(wt){
        compr_frq[[i]] <- list(fish=fishery,len=out.len,wei=out.wei,
                                         catch=ann.catch,effort=ann.effort,cpue=ann.cpue)
      } else {
        compr_frq[[i]] <- list(fish=fishery,len=out.len,catch=ann.catch,
                                  effort=ann.effort,cpue=ann.cpue)
      }
    }
  }
  return(compr_frq)
}
